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ABSTRACT 



A simple physical model of a copper crystal surface was de- 
veloped. Atoms were considered as quasi-hard spheres which 
occupied perfect lattice positions. A computer simulation, based 
on energy consideration only, using the Monte Carlo method was 
developed, tested and used to study equilibrium surface micro- 
states. As a result of this study, four conclusions were drawn: 

1. This model holds promise for further investigation of real 
crystal surface phenomenon. 

2. Minimum energy considerations cause atoms to align them- 
selves preferentially in a (lio) direction on the surface of a 
face-centered cubic crystal. 

3- Stepped surface configurations are fairly stable, but 
isolated "stub" atoms and vacancies tend to coalesce with other 
"stubs" and vacancies, respectively. 

4 . Random motions of the individual atoms cause aggregates of 
atoms to break apart and recombine. 
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I . INTRODUCTION 



A. HISTORY 

This simulation was based on two areas of contemporary research 
crystal growth and radiation damage. Ideas were borrowed from each 
of these fields of study in order to devise a new, simple physical 
model of a real crystal surface, and to begin the study of crystal 
surface equilibrium microstates. Brief discussions of prior re- 
search is contained in the following paragraphs. 

Theoretical models of crystal growth have been studied for over 
20 years, and numerous papers have been written attempting to de- 
velop models for growth from a solution, a melt, or a vapor 
[4,6,9,10] . An essential part of these investigations was the 

assumption of a perfect simple cubic crystal lattice and of first 
nearest neighbor forces only. Some results have been compared with 
experimental work done on the growth and dissolution of simple 
cubic ionic crystals in supersaturated and unsaturated solutions, 
and have shown both agreement and disagreement with the theory [6] . 
Motion of atoms absorbed on the surface of the crystal were taken 
into account by incorporation of theoretical surface diffusion mech 
anisms based on transport theory. These papers were generally 
written in terms of classical thermodynamic quantities. Crystal 
growth, diffusion and dissolution were all discussed in terms of 
Gibbs free energy and chemical potential for generalization to any 
simple cubic crystal. The crystals visualized were usually ionic. 
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Radiation damage studies provided useful physical information 
on real crystals. Gibson, Goland, Milgram and Vineyard (hereafter, 
GGMV) , and Girifalco and Weizer (hereafter, GW) have investigated 
interatomic potentials which were found useful. GGMV used a Born- 
Mayer exponential approximation to the interatomic potential, and 
applied their results to radiation damage dynamics for short in- 
teratomic spacings C 53 - The GGMV potential used here has become 
known as the "Gibson II potential". GW used a Morse potential and 
fitted data for 15 cubic metals to arrive at a appropriate Morse 
potential parameters for metals [ 7 ]. Anderman followed the GW 
procedure and obtained slightly different parameters for use in 
the Morse potential, by introducing the idea of truncating the 
potential function. Anderman used his values in some static radia- 
tion damage simulations [l]. The Anderman potential values have 
been used here, and the Anderman form of a GW Morse potential will 
be referred to as a "truncated Morse potential". 

Since solving for the motion of the atom in a crystal involves 
solving an N-body problem, an essential assumption for the radia- 
tion damage studies was that the motion of an atom could be des- 
cribed by an independent set of two-body potentials such as the 
Gibson II and truncated Morse potentials. Such a technique yielded 
a soluble approximation to the actual physical problem. For ex- 
ample, an atom and four of its first nearest neighbors gave four 
independent two-body problems, instead of one five-body problem. 

B. THE NEED FOR NUMERICAL SIMULATION 

The detailed processes involved in both of these areas of study 
were essentially unknown. Both lines of research were attempts to 



10 



take macroscopic information and develop microscopic models which 
predict macroscopic phenomenon. In some cases the theory involved 
did not permit analytical solution of the problem [6] 
and in all cases random processes at the atomic level were assumed. 
Solving the same problem a number of times using a Monte Carlo model 
and averaging the properties of several final configurations was a 
technique common to both of these areas of interest. Since the 
calculations were repetitive and could be programmed fairly simply, 
computer simulation was a logical tool to use. In addition, since 
microscopic models were assumed and system microstates were gene- 
rated, a computer simulation yielded microscopic information - 
which may be appropriately averaged and interpreted to yield theo- 
retical macroscopic results for comparison with empirical macro- 
scopic results. A "computer experiment" therefore yielded 
microscopic as well as macroscopic information. 

C. THE MONTE CARLO METHOD 

The Monte Carlo method requires that in calculation, whenever 
one makes a decision based on probability, the decision is made by 
generating a random number and comparing it with the probability 
involved. The calculation is continued until it arrives at another 
decision based on probability. The method then requires the gene- 
ration of a new random number for comparison and decision. The 
entire process continues for as long as desired, and a set of final 
conditions is obtained from a set of given initial conditions. The 
entire process may then be repeated, starting with the original set 
of initial conditions and new random numbers, to arrive at a new 
set of final conditions. In the end there is an ensemble of final 
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states generated from a single initial state, and it is possible 
to compute the average properties of the final state C 12 ] . The 
average final state properties may be compared with experimental 
data to yield information on the quality of the theoretical micro- 
scopic model. Essentially, the method attempts to generate macro- 
scopic information from a microscopic system model. 

D. THE MODEL ATTEMPTED HERE 

The present effort is believed to be the first attempt to study 
the surface dynamics of an actual metal crystal affected only by 
the temperature of the system. No attempt was made to include 
growth or dissolution (adding or subtracting atoms) of a crystal. 
Rather, a simple physical approach was tried in an effort to study 
only the surface of an actual metallic crystal. What was visualized 
was a section of a perfect copper crystal surface that was initially 
’’piled" in different ways. Observations were made of the subsequent 
surface configurations (microstates) resulting from Monte Carlo 
calculations. Interatomic potentials were borrowed from the radi- 
ation damage studies. Transition probabilities were borrowed from 
the thermodynamic crystal growth studies. An appropriate computer 
program was written to simulate the copper surface dynamics, and 
several computer runs were made to obtain a rough idea of what 
appeared to be equilibrium microstates. The programming language 
used throughout this study was FORTRAN IV, and the computer used 
was the IBM 360/67. 
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II. THE SIMULATION MODEL 



A. ASSUMPTIONS 

Since this was a first attempt at solving a new problem, 
several simplifying assumptions were made. Perhaps the greatest 
assumption made in this calculation was that metal crystal binding 
energies could be computed from a single composite potential energy 
function, with the interactions between atoms considered as in- 
dependent two-body problems. The essential feature of this as- 
sumption was to treat the atoms of a crystal lattice as quasi-hard 
spheres which interacted by some "unknown" mechanisms to yield 
observed crystal behavior. 

The only metal visualized was copper , which forms a face- 
centered cubic (fee) lattice. For simplicity, the (001) crystal 
orientation was assumed. Lattice sites were visualized as ordered 
real triples in the positive octant of a rectangular coordinate 
system. The coordinate planes and axes bordering the positive 
octant were considered as belonging to the positive octant. See 
Figure 1. Although a rectangular parallelepiped was always gene- 
rated to provide possible lattice locations (hereafter called the 
"active lattice volume"), appropriate changes in the lattice gene- 
rator insured that the initial configuration of the surface could 
take any desired shape with the sole precondition that the only 
sites available were perfect lattice sites. The active lattice 
volume was considered as merely a framework of sites that were 
available for atoms- whether the sites were actually occupied or not 
depended on the details of the raicrostate. 
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In addition, it was assumed that only the top several atomic 
layers of the surface provide interesting results - in effect, this 
meant that the active lattice sites generated were placed on top of 
a perfect lattice substructure which was not allowed to interact 
with the surface sites in any way other than to be used in calcu- 
lating surface site potential energy. This further implied that 
above the active lattice sites all locations were always kept 
vacant . 

Several additional assumptions were made to simplify the model. 
These were: 

1. Atomic vibrations and zero-point energy were neglected. 

2. Lattice relaxation due to vacancies was neglected. 

3. No interstitial atoms were allowed. 

4. No impurity atoms were allowed. 

5 . Only even numbers of planes in the principal coordinate 
directions were allowed. 

Further discussion of these assumptions is contained in Appendix A. 
B. THE COMPUTER PROGRAM 

A computer program that provided an adequate picture of micro- 
scopic surface activity had to be written and tested. Overall 
program flow is shown in Figure 2. The main features of the pro- 
gram are discussed in the remainder of this chapter. Further 
details are provided in Appendix A. 

1 . The Lattice Generator 

The first task involved generating an appropriate lattice- 
work of crystal sites in the computer. For the (001) orientation 
of a face-centered cubic crystal (hereafter, "fee, (001) crystal") , 



available sites are described by coordinates such that the sum of 
coordinates is either an even or an odd integer for all positions 
in the crystal. The even choice was used in this simulation. 

After a valid lattice site was located, it was given a 
number as a label. Then an occupation index was defined for the 
site to retain a knowledge of the surface configuration. Several 
different initial configurations were used, and these are sketched 
as the "a M parts of Figures 3 through 17* The specific program 
statements to generate such surfaces are contained in Appendix B. 

2 • Periodic Boundary Conditions 

Computer storage and time limitations made impractical the 
use of an arbitrarily large active lattice volume. A practical 
lattice size involved at most a few thousand active lattice sites. 
This required that the active volume size be limited to a 
Cartesian space of roughly 20 x 20 x 10, which corresponds to 
2,000 active lattice sites and a copper crystal measuring about 
5oX x 508 x 258 . Such a small crystal with six surfaces compli- 
cated the study of surface effects. 

One surface (the bottom of the active volume) has already 
been eliminated by assuming that the active lattice volume was 
superimposed on a perfect crystal substrate which was never allowed 
to vary. Four other surfaces (the sides of the active volume) were 
eliminated by assuming ’’per iodic boundary conditions” which resulted 
in an essentially infinite single surface (the ’’top” of the active 
volume) that varied periodically in configuration. This infinite 
surface may be visualized as part of a crystalline surface far 
from a grain boundary. 
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These periodic boundary conditions were achieved by 
visualizing that each active volume plane was laterally surrounded 
by exactly similar active volume planes that had been translated in 
the x and/or y dir ect ion( s ) . In other words, an atom on the left 
edge of an active volume plane "sees" the right edge of that same 
active volume plane; an atom on the top edge "sees" the bottom edge 
of the same plane, and similarly for the other edges. 

Periodic boundary conditions were achieved in practice by 
defining nearest neighbor (NN) arrays for each site. These arrays 
were calculated by visualizing the spatial relationships between 
an atom in the fee, (001) crystal and its nearby surrounding sites - 
see Figures 18 and 19. Atoms within one atomic layer of the sur- 
face of the active lattice volume required special attention to 
insure the proper periodic structure of each plane. 

3 • Potential Energy 

The model chosen assumed that the crystal may be approxi- 
mated by a latticework of quasi-hard spheres, so a scheme for 
holding individual atoms together was needed. The problem was an 
N-body problem, and therefore no general exact solution existed. 
However, approximating the potentials by an independent set of two- 
body potentials yielded a solvable approximation to the physical 
system, and this was a standard technique used in radiation damage 
studies . 

The next task involved a judicious choice of potential 
function. A composite potential function was chosen. For low 
values of interatomic spacing the Gibson II potential C 53 was used, 
and for larger values the truncated Morse potential Cl] was used. 
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Since these two potentials do not match smoothly together (see 
Figure 20), a cubic function was generated which matched the value 
and the slope of these two potentials at appropriate points so 
that a smooth composite potential was achieved (see Figure 21). 

The cubic potential coefficients were generated by solving a system 
of four simultaneous linear equations using Subroutine CROSYM 
developed at Lockheed. The form of the potential equations, their 
appropriate coefficients, and the cutoff values used in this simu- 
lation are listed in Table I. 

Next, a potential energy was associated with each site in 
the lattice based on the number of its first and second nearest 
neighbor (NNl and NN2, respectively) sites that were occupied. In 
other words, every active volume lattice site was given a potential 
energy. If the site were occupied, this energy was just the total 
mutual interaction energy between an atom and its NNl's and NN2 ' s . 
However, if the site were vacant, this energy was the interaction 
energy between the NNl's and NN2's of the site that would arise if_ 
the site were occupied. In practice, potentials were given to each 
site (occupied or vacant) in the same manner , but the use to which 
these potentials were put depended vitally on whether a particular 
site was occupied or not. 

On the average about half of the mutual interaction energy 
between sites will be associated with each site individually, so a 
parameter PFIV (mnemonic for "point five", and usually equal to 
0.5) was used as an energy distribution factor. By varying this 
energy distribution factor, physical effects resulted which were 
equivalent to varying the lattice temperature. Therefore, by 
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changing the energy distribution factor, different lattice tem- 
peratures were achieved. The relationship between the equivalent 
lattice temperature and PFIV is indicated in Figure 22 and dis- 
cussed in greater detail in Appendix A. 

4 . Transition Probability 

If a given lattice site were occupied and had one or more 
vacant sites as NNl's it had (statistically) a finite probability 
of moving from its present site to the nearby vacant site. In 
practice, when an occupied lattice site was reached, the computer 
searched for lower values of site potential among the unoccupied 
NNl 1 s of the occupied site of interest. Such lower energy vacant 
sites implied that the atom under consideration was sitting on the 
edge of a potential well, and a transition was carried out to the 
deepest NNl well - that is, the atom " jumped” into the deepest 
NNl well. If a "deepest well” was not unique, then the actual 
final site was chosen stochastically. "Wells” of zero depth were 
considered as having "lower energy" for the purpose of this 
calculation . 

For atoms with no vacant NNl sites of lower energy but 
one or more with higher energy, the atom would have to climb a 
potential hill to reach it (them). This was considered possible 
but was given a probability defined by a Boltzmann factor: 

P = exp (-E/kT) (1) 

where 

P = the probability of transition to a higher energy site 
E( > 0) = energy difference between the occupied site of 

interest and one of its NNl sites with higher energy 
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k = Boltzmann’s constant 



T = absolute temperature 

The computer searched for the lowest positive energy difference 
and made the decision to jump or not stochastically. If several 
NN1 sites had positive energy differences, they were considered 
in turn from the lowest to the highest until a transition was made. 
If no transition were made by the time the positive energy dif- 
ferences were exhausted, then the atom stayed where it was for 
that particular microstate. 

After any transition the potential energy of the NNl's and 
NN2*s of both the old and new sites were adjusted to account for 
the changes made. After every third microstate, all potentials 
were recalculated to avoid loss of accuracy. A pseudo-random 
number generator was used to produce the random probabilities 
needed to make the decisions involved. 

5. Computer Output 

The problem of how to present microstate information in 
easily interpreted form was solved by having the computer ’’draw” 
pictures of microcrystal planes, called ’’microstate pictures”. 

The computer printed arrays of occupation indices to indicate a 
vacancy or an occupied site, and each digit was in its correct x 
and y coordinate position for the site it represented. See 
Figure 23. 
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III. RESULTS 



A. INITIAL TESTING 

Test runs made prior to taking any ’’good data” indicated that 
two problems would be encountered. The first problem discovered 
was the lack of net motion out of a perfect surface configuration. 
The second problem was an indication that the model might tend to 
"transport" atoms preferentially in any given plane from higher to 
lower values of the arbitrarily assigned y-coor dina te. 

The problem of no net motion from a perfect surface was solved 
by two different techniques. One was to alter the energy distri- 
bution factor PFIV to define an equivalent lattice temperature to 
literally force some action. The other method involved piling the 
surface atoms in both random and preconceived configurations for 
the initial microstate and observing subsequent microstates using 
the correct average value of PFIV (= 0.5). The results of these 
tests are discussed in Sections B and C, below. 

The problem of model dependent transport was more fundamental, 
since any results obtained would be biased in some unknown fashion 
Again, two techniques were used to explore the problem. The first 
technique was to test the motion of a single atom initially placed 
on a perfect surface. The motion of this single "stub" atom was 
observed for 1+00 raicr ostates . The motion within each microstate 
was recorded as the net number of lattice units moved "up", "down" 
"left" or "right" from its previous position. These directions 
were defined for an observer looking at a single lattice plane; 
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"right" meant "in the + x-dir ect ion " , and "up" meant "in the 
+ y-direction" for the plane under consideration. The average 
results are reported in Table II to the nearest one percent. A 
quick glance at the figures appeared to lead to the conclusion that 
the model was biased to motion right and down. However, consider- 
ation of the average number of jumps weighted by the distance 
jumped indicated that: 

1. When the stub moved left, its average motion was 0.47 
lattice units . 

2. When the stub moved right, its average motion was 0.6l 
lattice units . 

3. When the stub moved up, its average motion was 0.46 
lattice units . 

4. When the stub moved down, its average motion was 0.44 
lattice units . 

This meant that the average net motion of a single stub atom was 
right and almost entirely lateral. 

Since this test was performed with only one stub atom, a second 
technique was used to test for model biasing. The technique was 
to study the apparent motion, if any, of large numbers of atoms 
placed on a perfect surface. This required, in effect, "taking 
data" and the results were available only after several long com- 
puter runs were made. Based on the runs discussed below, it was 
concluded that there existed no significant motion of large numbers 
of atoms attributable to model biasing. Some motion "down" was 
found, but this amounted to about one lattice unit every 10 to 20 
microstates for 100 or more atoms and was not considered 
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significant. This observation was purely qualitative, and such 
motion was not found in every case. In addition, no gross lateral 
motion was observed, and any that did exist was much slower than 
the "down" motion. 

B. PERFECT SURFACES 

Tests made using perfect surfaces required that the energy 

distribution factor PFIV be altered. The arbitrary initial tem- 

o 

oerature chosen for all calculations was 1,000 K. With the 
correct average value of PFIV (= 0.5), no net change was observed 
for any microstate in 500 attempted microstates when the initial 
microstate was a perfect surface. 

In an effort to force some different microstates to occur, 

PFIV was lowered by a factor of 10, yielding an equivalent tem- 
perature of 10,000 °K - well above the copper boiling point of 
2,583 °K . This resulted in a decidedly unphysical sequence of 
surface microstates. Since the temperature was so high, the sur- 
face literally evaporated and burroughed a deep hole into itself 
within 10 microstates. The model restraints insured that no atom 
would be lost, but this did not mean that the microcrystal did not 
try to "explode" - a program flag included to warn that an atom had 
no NNl f s was flashed at least 35 times during simulation of 500 
microstates and never appeared in any other simulation. Equilibrium 
numbers of atoms in each plane was reached after about 30 micro- 
states (as opposed to more than 100 microstates normally required 
for other simulations). Some single atom vacancies and stubs were 
observed, but such surface defects were repaired quickly. However, 
the frequency of observation of these defects was much greater in 
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this simulation than for all other simulations. Figure 24 shows 
a typical microstate achieved. Note that the figure is drawn with 
lines bordering the atoms and atom clusters in a (lOO) direction 
(Figure 24a) and in a (llO) direction (Figure 24b). For all other 
simulations, lines parallel to a (llO) direction appeared to be 
most natural to a face-centered cubic crystal, but here the jagged 
edges appearing in both sketches indicated that neither was 
"natural". Chapter IV contains a discussion of the "natural" 
crystal orientation. Also note that there are several "overhangs" 
visible in each sketch - overhangs were observed in other simu- 
lations, but not with the great frequency with which they appeared 
her e . 

The next perfect surface attempted was at an equivalent tem- 
perature of 1,355 °K - just one degree Kelvin below the copper 
melting point. In 2G0 microstates no net change was observed for 
any raicrostate. 

Next, a perfect surface at an equivalent temperature of 
2,000 °K was attempted. A typical microstate achieved is indicated 
in Figure 25. The atoms making up the top perfect surface layer 
"jumped up" and coalesced on the surface. 

Finally, a simulation at 2,584 °K (one degree Kelvin above the 
copper boiling point) was attempted. Equilibrium numbers of atoms 
in each plane were reached by about microstate 100. A typical 
microstate is indicated in Figure 26. Atoms in several surface 
layers jumped up and actually formed steps with edges parallel to 
the (llO) direct ion o 
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C. PERFECT SURFACES PLUS EXTRA ATOMS 



Most simulations attempted fell into this category. This type 
of initial configuration was used for two reasons: 

1. It was unphysical to expect that real crystal surfaces 
would be absolutely perfect for large sections of surface. 

2. It was unphysical to require temperatures greater than the 
melting point in order to achieve some idea of what surface equili- 
brium raicrostates were. 

All of these simulations were run at a lattice temperature of 
1,000 °K and using the correct average value of the energy distri- 
bution, PFIV (= 0.5) . 

Two simulations were run placing atoms randomly on top of a 
perfect surface. Different random number seeds were used for each 
simulation, resulting in different numbers of extra atoms on top 
of a perfect crystal surface. See Figures 3a and 4a. The randomly 
placed atoms coalesced within 20 microstates in both simulations. 
Figures 3b and 4b show typical equilibrium configurations for each 
experiment. Note that the atoms arranged themselves preferentially 
along the (llO) direction. 

One simulation was run placing an extra half plane of atoms on 
top of a perfect surface, with the halfplane edge parallel to the 
x-axis. See Figure 5a. At equilibrium, the atoms were still to- 
gether but again showed a preference for lining up parallel to the 
(llO) direction. See Figure 5b. 

Two simulations were run piling atoms into a pyramid config- 
uration with sides of different slopes. Sketches of the pyramid 
initial microstates may be found as Figures 6a and 7a with the 



24 



atoms aligned in the ( 100) and (010) directions. The final con- 
figurations are shown in Figures 6b and 7b. Note that in Figure 6b 
all piled up atoms have fallen down to coalesce on top of the 
original base plane. Note that in Figure 7b not all of the atoms 
have fallen down - in fact, there are aggregates of atoms in four 
planes above the original base plane. The large "hole" appearing 
in the plane above the base plane is the outer "groove" appearing 
in the original microstate that has been reoriented and partially 
filled in. Note that one atom has actually jumped out of the base 
plane and has not been replaced - this is represented by the "hole" 
in the left corner of Figure 7b. 

Two simulations were run with the atoms arranged to form a 
"pyramid hole" of different step widths - see Figures 8a and 9a. 

The configurations 500 raicrostates later are shown in Figures 8b 
and 9b. Note that in both cases the atoms have fallen down into 
lower planes, but that the "falling down" process is a bit slower 
when the hole step is 2 atomic layers wide then it is for a hole 
step of 3 atomic layers. The reason for this is probably the fact 
that with 2 atomic layers per step there are more atoms in the 
active lattice volume, and the active volume itself is larger. 

Three simulations were run piling atoms into ridges - one was 
parallel to the x-axis and the other two were parallel to the 
y-axis. See Figures 10a and 11a. Despite the different ways the 
atoms were initially piled and despite the fact that different 
random number seeds were used for the two ridges paralleling the 
y-axis, there was a remarkable resemblance in the gross appearance 
of equilibrium microstates of each experiment - these are sketched 
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as Figures 10b, lib and 11c. The slots appearing in these sketches 
are all in the x = 1 plane, which was not completely filled in the 
initial raicr ostates , i.e., the holes do not represent several atoms 
jumping out of the z = 1 plane but the rearrangement of a hole 
originally placed there after it has received some additional atoms 
from the higher planes. 

Next, four simulations were run in which the atoms were piled 
into valleys of different step widths. The initial microstates 
are sketched in the "a" parts of Figures 12 through 15. Note the 
similarity between Figures 10a and 12a and between Figures 11a and 
l4a “ in effect, since periodic boundary conditions were applied, 
Figure 12a may be obtained from Figure 10a by translating Figure 10a 
one-half the active volume distance in the y-direction, and simi- 
larly for Figures 11a and 14a with translation in the x-direction. 
In any event, it may be expected that some overall gross similarity 
of the final microstates may appear. Comparison of Figures 10b and 
12b, and Figures lib and l4b tends to verify this expectation. 

Note that again the holes appearing in the final microstates were 
a result of rearrangements and partial fill-in of the original 
holes in the 2=1 planes. Figures 13a and 15a indicate valleys of 
narrower step width, and final microstates achieved are sketched 
as Figures 13b and 15b. 

Since there was a preference for alignment along the (lio) di- 
rection, the two final simulations started with all atoms aligned 
in that direction. One was visualized as a truncated (1 -1 1) 
plane (Figure l6a), and the other as a truncated (-1 1 1) plane 
(Figure 17a) . The 500th microstates are sketched as Figures l6b 
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and 17b. Note that in each case atoms were still piled high and 
no plane was completely filled - these facts indicate that the 
system probably had not reached equilibrium. These initial micro- 
states were a bit more unphysical than the others tested since they 
started with steps several atomic layers high on 2 sides and one 
atomic layer high on the other . This permitted overhangs when 
atoms "jumped off" the steep sides, but this was not considered 
seriously detrimental to the results. 



27 



IV. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

As indicated in the discussion of the perfect surface simu- 
lation results, most "exper iments M yielded results which imply 
that the (110) orientation was most natural for sketching equilib- 
rium microstates. Here, "natural” meant that there were fewer 
jagged edges and protruding atoms than if the microstate were 
sketched in the (100) orientation. Note that the initial con- 
figuration sketches (the "a" part of Figures 3 through 17) are all 
in the (100) orientation, since this was how the crystal was initi- 
ally visualized and developed. The following discussion will 
neglect the one unphysical case achieved - that of an initially 
perfect surface at 10,000 °K. 

The first important result of testing this model was that there 
was one thread of similarity running through all the experiments - 
in the final microstates achieved, all atoms coalesced into steps 
either parallel or perpendicular to a (llO) direction. This was 
most encouraging since the initial state had no effect on the final 
step direction. It was further encouraging since when building an 
fee, (001) crystal from marbles and glue, atoms on the surface are 
more closely arranged in a (llO) direction than in a (lOO) direction. 
The physical reason for such a preference lies in energy considera- 
tions: in an fee, (001) crystal, NN1 f s are located along (lio) 
directions, and NN2*s along (lOO) directions. Therefore, the 
lower energy state sought by a system intuitively and deliberately 
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built into the transition probabilities of this model indicate a 
preference for the (llO) direction. Again, it is a pure energy 
consideration which dictated that the atoms coalesce and not dis- 
perse randomly. The minimum energy condition was achieved for as 
many atoms as possible staying close together. Occasionally, iso- 
lated surface vacancies or stub atoms would appear, but they were 
quickly repaired or moved so that these defects tended to coalesce 

The second conclusion involved the final shape of piled atoms. 
Atom piles tended to be stable in a stepped configuration over 
several microstates, and this agreed with preconceived notions of 
surface shape. Figure 27 shows a sketch due to D.E. Harrison [8] 
indicating his concept of a surface equilibrium microstate. Ac^ 
cording to the model developed and tested here, the sketch is 
essentially correct if it is drawn in a (110) orientation for its 
stepped features. However, this model indicated that isolated 
holes and stubs do not persist at equilibrium. Isolated vacancies 
or "holes" were as mobile as isolated stub atoms. Isolated groups 
or atoms and holes appeared to move slowly as a group to coalesce 
with larger aggregates. 

The third observation of importance involved breaks and their 
repair in elongated strings of surface atoms. Figure 28 shows a 
series of sketches illustrating this phenomenon. When the motion 
of the atoms led to the movement of atoms in opposite directions, 
a local '’necking down” of surface aggregates of atoms was observed 
Random motions led to selective breaks at such thin steps, but 
unless the atoms in the two resulting steps moved several atomic 
diameters apart the break actually repaired itself. Conversely, 



29 



unless the break repaired itself by supplying several atomic 
diameters of step, the break was likely to recur within 10 or 
20 raicr ostates . This behavior was observed repeatedly, and was 
not really surprising although it was not expected beforehand. 

B. RECOMMENDATIONS 

First, it is recommended that this model be further studied 
and refined to obtain a more realistic picture of metal surface 
microstates. The results achieved and conclusions drawn here must 
be considered only as preliminary indications since this was a 
first attempt at a new physical model and since the testing was 
not extensive. 

Secondly, it is recommended that prior to extensive experi- 
mentation to find equilibrium microstates, the problem of model 
dependent transport be more thoroughly investigated. Since large 
aggregates of atoms provided the most interesting results here, 
this testing should be carried out with groups of at least 10 
atoms on a perfect surface. It is envisioned that a practical test 
might involve calculating the center of mass of these aggregates 
for several hundred microstates and computing the average motion 
of their center of mass. If trends appear (as in the case of a 
single stub atom, tested here), then appropriate corrections must 
be made in the computer program. Perhaps the easiest possible 
change to install involves the sequence in which the atoms are 
looked at as a microstate is generated. As it was used here, 
every microstate was generated by having the computer investigate 
each atom in order by lattice site number, which meant that r ows 
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were investigated from left to right and from bottom to top in 
each z-plane. It may be more realistic to generate every other 
roicrostate in a different way, e.g., by investigating each column 
from bottom to top and from left to right in each z-plane. 

Third, it is recommended that the following changes be made in 
order to improve the model: 

a. Rewrite this fee program to generate the active lattice 
volume in a (110) orientation. This would require a new lattice 
generator, new lattice loading steps, new NN assignment segments, 
and recalculation of potential energies. 

b. Include consideration of NN3 t s and NN4 1 s - this would 
require new input parameters for the truncated Morse potential 
plus NN3 and NN4 assignment segments. 

c. Allow floating point values for the atom positions and 
introduce kinetic energy for the atoms. It is expected that 
introducing a normal distribution of kinetic energies would re- 
duce the total energy of perfect surface atoms sufficiently to 
start them "jumping" out of their positions without resorting to 
the equivalent temperature idea introduced here. This would al- 
low for the introduction of lattice temperature as an input para- 
meter and controlling factor in lattice activity. 

d. Write a lattice loader that would generate a more realistic 
lattice - e.g., one in which the atoms were randomly placed sub- 
ject to physical restraints such as having at least one NN1 and 
occupying nearly perfect lattice positions (i.e., this assumes 

that floating point lattice positions have already been 
incorporated) . 
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e. Allow for vacancies, interstitials and impurities. This 
would require consideration of lattice relaxation and strain 
effects, and would require new potentials for impurity-copper and 
impurity-impurity interactions. 



32 



APPENDIX A 



THE DETAILS OF THE SIMULATION MODEL 
A. ASSUMPTIONS 

Perhaps the greatest assumption made in this calculation is 
that metal crystal binding energies can be computed from a single 
composite potential energy function, with the interactions between 
atoms considered as independent two-body problems. The essential 
feature of this assumption was to treat the atoms of a crystal lat- 
tice as quasi-hard spheres which interact by some "unknown" mecha- 
nisms to yield observed crystal behavior. The "unknown" mechanisms 
include, in the general case, attraction due to ionic, covalent, 
metallic and van der Waals forces, and repulsion due primarily to 
electron cloud overlap. The specific mechanisms, however, were of 
no concern. No quantum mechanical considerations were attempted - 
that is, the calculations involved were all classical once the ap- 
propriate composite potential function was assumed. First and 
second nearest neighbor interactions were considered. 

The only metal visualized was copper , which forms a face- 
centered cubic (fee) lattice. For simplicity, the (001) crystal 
orientation was assumed. Lattice sites were visualized as ordered 
real triples in the positive octant of a rectangular coordinate 
system. The coordinate planes and axes bordering the positive 
octant were considered as belonging to the positive octant. See 
Figure 1. Although a rectangular parallelepiped was always gene- 
rated to provide possible lattice locations (hereafter called the 
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"active lattice volume") , appropriate changes in the lattice gene- 
rator insured that the initial configuration of the surface could 
take any desired shape with the sole precondition that the only 
sites available were perfect lattice sites. The active lattice 
volume was considered as merely a framework of sites that were 
available for atoms - whether the sites were actually occupied or 
not depended on the details of the microstate. 

It was further assumed that atomic vibrations about lattice 
sites could be neglected, so there were no zero-point energy con- 
siderations or fluctuations in lattice potential due to atomic 
motion. Lattice relaxation due to crystal vacancies was ignored, 
as were interstitial atoms and their consequent lattice strain 
effects. Furthermore, no impurity atoms were allowed. With these 
assumptions, it was possible to describe each available lattice 
site as ordered integer triples and use a single composite po- 
tential energy function for every two-body pair in the system. 

This simplified the model considerably and yielded a great re- 
duction in the potential energy calculations required. 

In addition, it was assumed that only the top several atomic 
layers of the surface would provide interesting results - in 
effect, this meant that the active lattice sites generated were 
placed on top of a perfect lattice substructure which was not 
allowed to interact with the surface sites in any way other than 
to be used in calculating surface site potential energy; this fur- 
ther implied that above the active lattice sites all locations were 
always kept vacant. Specifically, latt-ice sites described by z<0 
were considered always full - see Figure 1. Lattice sites 
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described by O ^ 2 ^ z were considered as available for con- 
max 

taining an atom, depending on the particular microstate initialized 

or achieved during a calculation. Lattice sites described by 

z > z were considered always empty, 
max 

Finally, this model was restricted to even numbers of planes 
in each of the principal coordinate directions to simplify the 
calculations involved in assigning first and second nearest neigh- 
bors to each available site. This means that x , y , and z 

max max 7 max 

must all be odd, since each coordinate plane was included as part 
of the active lattice. 

B. THE COMPUTER PROGRAM 

1 . The Lattice Generator 

a). Logic. Before any calculations could be made, the com- 
puter had to be told what the available sites for atoms were and 
whether or not these sites had an atom - that is, whether a given 
site was occupied or vacant. For the ( 001 ) orientation of an fee 
crystal (hereafter, M fcc, (001) crystal"), the essential technique 
involved relied on the fact that the sum of the site coordinates 
must be an even (or an odd) integer (this point is particular to 
this lattice shape and orientation). Refer to Figure 29. For 
example, the Cartesian point (1, 1, 1) has a sum of coordinates 
equal to three, an odd integer - therefore, this point is not an 
acceptable lattice location. On the other hand, Cartesian points 
(1, 0, 1) and (1, 1, 0) each have a sum of coordinates equal to 
two, an even integer - these points are acceptable lattice locations 
and are assigned integer labels when they are reached in the lattice 
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generator segment of the program. Figure 30 shows the numbering 
scheme for lattice points. 

The lattice was generated by a set of nested DO loops. 
One plane at a time was generated, and these planes were parallel 
to the x-y plane. The coordinates of each acceptable site and its 
occupation index were stored as vector arrays. An ’’occupation 
index” (vector ”N0CC(I)”) merely indicated whether or not a site 
was occupied: an index of zero implied that the site was vacant, 
while an index of one meant that it was occupied. No values other 
than zero or one were allowed for occupation indices.' 

b) . The Initial Microstate. In order to provide adequate 
testing of the simulation model, the initial microstate must be 
easily changed. This was achieved by substituting program segments 
containing one or more cards into the program. Initial microstates 
actually used were: 



1) 


Perfect surface 




2) 


Perfect surface plus extra atoms 






( a ) 


randomly dispersed extra atoms 


- Figures 






3a and 4a 






(b) 


a monatomic step - Figure 5a 






(c) 


pyramids - Figures 6a and 7a 






(d) 


pyramid holes - Figures 8a and 


9a 




(e) 


ridges - Figures 10a and 11a 






(f) 


valleys - Figures 12a, 13a, l4a 


and 15a 




(g) 


(1 ~1 1) plane - Figure l6a 






(h) 


(-111) plane - Figure 17a 





The actual program segments which generate these different initial 
conditions are listed in Appendix B. 
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A site labeled "IFUL" was used to represent all sites 

in the crystal substrate, that is, all sites for z < 0; since these 

sites were always full (occupied), NOCC(IFUL) was always one. A 

site labeled "IVAC" was used to represent all sites above the 

active lattice volume, that is, for z > z ; since these sites 

’ max 

were always vacant, NOCC(IVAC) was always zero. These assignments 
were made immediately after the lattice generator and were based 
on the last value of ,r M" reached in the lattice generator. 

2 . Periodic Boundary Conditions 

a) . The need. Computer storage and time limitations made 
impractical the use of an arbitrarily large active lattice volume. 

A practical lattice size involved at most a few thousand active 
lattice sites. This required that the active volume size be limit- 
ed to a Cartesian space of roughly 20 x 20 x 10, which corresponds 
to 2,000 active lattice sites and a copper crystal measuring about 
50 X x 50 X x 25 X. Such a small crystal with six surfaces com- 
plicated the study of surface effects. 

b) . The method. One surface (the bottom of the active 
lattice volume) was eliminated by assuming that the active lattice 
volume was superimposed on a perfect crystal substrate which was 
never allowed to vary. Four other surfaces (the sides of the 
active volume) were eliminated by assuming periodic boundary con- 
ditions which resulted in an essentially infinite single surface 
(the 1T top" of the active volume) that varied periodically in con- 
figuration. This infinite surface may be visualized as part of a 
crystalline surface far from a grain boundary. 
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The overall scheme was represented as follows (see 
Figure 31) : Each active volume plane was laterally surrounded by 
exactly similar active volume planes that had been translated in 
the x and/or y dir ection(s ) . In Figure 31a, the actual active 
volume was visualized as the center square; the remaining eight 
squares were the "boundary volume planes" which had exactly the 
same configuration as the center square. Square 2 was achieved 
by translating square 1 a full microcrystal distance in the "-x" 
direction, and square 3 was achieved by translating square 1 a 
full microcrystal distance in the "+x ft direction. Squares 4-5-6 
were achieved by the similar process of translating squares 
2-1-3 a full microcrystal distance in the f, +y" direction, and 
similarly for squares 7-8-9. 

The active volume z-plane of Figure 31b is por- 
trayed as having six x-planes and six y-planes, and a total of l8 
lattice sites; three of these lattice sites are labeled (A,B,C), 
and the location of the remaining sites are indicated by X's. The 
images of the labeled sites are indicated in each of the surrounding 
boundary planes. Without the boundary planes, a vector [l,l,o] 
from active volume site "C" would yield a point in free space. 
However, with the boundary planes, as in Figure 31b, that same 
vector from "C" yields a site labeled "B" in the boundary plane. 

In effect, looking to the right from the active volume plane, one 
passes through infinity and comes upon the left edge of the same 
active volume plane in the same row from which one started. Simi- 
larly, a vector [o,6,o] from site "A" brings one right back to 
site "A". (This is merely an illustration, and actual crystals 
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generated typically had 20 or so planes in the x and y directions, 
so that returning to the same site was not so easy in practice.) 
During the calculation, any image sites reached were treated as 
active volume sites. Specifically, if in the active volume, sites 
A and C were occupied, but site B was vacant, then in all eight 
boundary planes, sites A and C were occupied and site B was vacant. 
This gave a periodic character for site occupation as one moved 
across a given z -plane, hence the name periodic boundary conditions. 
In actuality only one active volume was used, but the "image site” 
and "boundary plane" ideas explained here help to illustrate how 
periodicity was achieved. 

c) . Second nearest neighbor assignment segment. The NN2 
assignment segment is easier to understand than the NN1 assign- 
ment segment, and so is considered first. Figure 19 illustrates 
the positions of the six NN2 ! s in a fee, (001) crystal, using the 
numbering scheme developed for the program. To find the site 
number of any second nearest neighbor of a given atom, the given 
atom was visualized as the center of coordinates of a small mobile 
coordinate system as depicted in Figure 19. In general, the site 
at the origin was labeled "I", and the six positions for NN2 T s were 
labeled "J", where J varied from one to six. A double subscript 
notation was chosen and named array "NBRTW3 ( J , I) " - a mnemonic for 
"neighbor two" - where subscript J refers to one of the six NN2 
positions, subscript I refers to the label of the lattice site 
being investigated (given to that site by the lattice generator), 
and the value stored at the particular position was the label of 
the lattice site at position J with respect to site I. 
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In the general case, NN2 position one was labeled 
site number (1-1), and NN2 position two was labeled site number 
(1+1). NN2 positions three and four required that a shift of two 
rows of atoms in the "-y" and "+y" directions be made. For even 
numbers of x planes, that is, for IX even, the appropriate number 
of site numbers to "shift” is IX, since each row of sites in a 
given plane contains IX/2 sites (IX even). Consequently, NN2 
position three was labeled site number (I-IX) and NN2 position four 
was labeled site number (I+IX) . NN2 positions five and six re- 
quired a shift of two planes of atoms in the "-z" and "+z" di- 
rections. Since each row of a given plane contained IX/2 sites 
(IX even), and there were IY such rows per plane, each plane con- 
tained (IX/2)*(IY) sites. Therefore, two such planes contain 
IX*IY sites. As a result, NN2 position five was labeled site 
number (I-IX*IY) and NN2 position six was labeled site number 
( I+IX*IY) . 

For sites on, or one atomic layer in from, the active 
volume surface, proper labeling of NN2 sites helped achieve the 
desired effect for a single infinite surface varying periodically. 
For example, a site with NX(I)=0 (on the left face of the active 
volume) or NX(I)=1 (one atomic layer in from the left face of the 
active volume) looking for its neighbor at NN2 position one must 
"see" that neighbor (call it I’) with NX(I ? )=z -1 or NX(lM=z , 
and NY(I’)=NY(I), and NZ ( 1 1 ) =NZ( I ) . In the site labeling scheme 
chosen, this represents a shift of (IX/2)-l site numbers (IX even). 
Similar shifts must be made when "looking through infinity" from 
the other three sides of the active lattice volume, and this was 
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essentially how the imposed periodic boundary conditions were 
achieved in practice. 

For sites with NZ(I)=0 or NZ(I)=1 looking for their 
neighbor at NN2 position five, site IFUL (the always full sub- 
strate position, was reached. Similarly, for sites with NZ(I)=z 
or NZ( I ) =2 max “ 1 looking for their neighbor at NN2 position six, 
site IVAC (the always empty above lattice position) was reached. 

In this way a perfect substrate below the active volume and a void 
above the active lattice were maintained. Further specific de- 
tails of the NN2 assignment segment are omitted here, but may be 
found by studying the actual program between statements number 
195 and 201. 

d) . First nearest neighbor assignment segment. The NNl 
assignment segment was approached in much the same way as the NN2 
assignment segment. The problem was complicated by the fact that 
there are 12 NNl’s for an fee, (001) crystal, and there are 22 
special cases (6 faces, 12 edges, and 4 corners) to consider, as 
opposed to six NN2 v s and six special cases for the NN2 assignment 
segment. Figure 18 illustrates the positions of the 12 NNl's, 
using the numbering scheme developed for the progrant. Again, a 
double subscripted array was chosen to store NNl information, and 
it was called T, NBR0NE( J , I) " - a mnemonic for ’’neighbor one” - where 
subscript J refers to one of the 12 NNl positions, subscript I re- 
fers to the lattice site under consideration, and the value at the 
particular array position is the site number of the site at position 
J with respect to site I. 



The initial complication that arose was an 1f even-odd" 
dependence on the x-value of the site for 8 out of 12 NN1 positions. 
Taking this into account, the overall problem was solved by first 
assigning each site its 12 NNl 1 s without regard to the 22 special 
cases. These were computed using a shift technique similar to 
that developed in the NN2 assignment segment. The 22 special cases 
involved sites on the faces, edges, and corners of the active lat- 
tice volume, and required particular attention in order to appro- 
priately achieve the desired periodic boundary conditions. The 
method used then involved testing for each of the 22 special cases. 
As each special case was found a series of from four to nine ad- 
justments were made from the "general NNl sites” already stored. 
Sites IFUL and IVAC were employed in a manner similar to that 
involved in the NN2 assignment segment. Further details of the 
calculation are tedious and will not be examined here. The actual 
procedure may be gleaned from the program itself, beginning after 
the lattice generator and running through statement 195- 

Thus, the periodic boundary conditions imposed on the 
active lattice volume were achieved in practice via the use of the 
NNl and NN2 assignment segments. These segments, as written, re- 
quired that IX, IY and IZ all be even. This was not considered a 
very severe restriction, and there is no a prior reason to expect 
that the results were biased due to a choice of even numbers of 
planes in the coordinate directions. 

3. Potential Energy 

a). The Composite Potential Energy Function. The model 
chosen assumed that the crystal may be approximated by a 



42 



latticework of quasi-hard spheres, so a scheme for holding indi- 
vidual atoms together was needed. The problem is essentially an 
N-body problem, and therefore no general exact solution exists. 
However, approximating the potentials by an independent set of 
two-body potentials yields a solvable approximation to the physical 
system, and was a standard technique used in radiation damage 
studies . 

Since the approach taken involved the solution of 
independent two-body potentials, the next task involved a judicious 
choice of potential function. The choice began with the Morse 
Potential, proposed in 1929 as an appropriate potential for the 
wave equation of a diatomic molecule C 13 3 r 

(A-l) V Morse = D ex P( _2a ( r - r 0 ^ “ 2D ex P< _a ( r ' r Q ) ) . 

where 

D = dissociation energy of the molecular bond 
r Q = equilibrium interatomic separation 
r = interatomic separation 

a = constant related to the curvature about the 
potential minimum of the function 

Note that the first terra dominates for small "r " , and that the 
second term dominates for large "r". 

Girifalco and Weizer calculated Morse potential con- 
stants for 15 cubic metals, and these have found satisfactory 
agreement with experiment C 73 • Anderman used essentially the same 
procedure as Girifalco and Weizer, but used a truncated potential 
and obtained parameters for a copper Morse-type potential for two, 
three, and four nearest neighbors Cl]. This "truncated Morse" 
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potential with parameters for two nearest neighbors was used for 

the range r < r < r , where r is an inner cutoff value chosen 
B C B 

as a lower bound for the truncated Morse potential, and r^ is the 
truncation distance. 

GGMV arrived at reasonable short range results using a 
Born-Mayer type exponential approximation for the interatomic 
potential : 

( A » 2 ) v G i bson II = D ex P(’ a ( r ~ r 0 )) ’ 

where D, a, r and r have the same meanings as in V 

’ * o v Morse 

This "Gibson II” potential was found to be a good approximation 

for small interatomic separations, and so was used for short 

distances, 0 < r < r, where r is an outer cutoff value chosen as 
A A 

an upper limit for the region of Gibson II potential validity. In 
general, r^ ^ r^ because these two potentials do not match smoothly 
- see Figure 20. The "mismatch region" was fitted by a series ex- 
pansion with a cubic function, since only four pieces of information 
were available to match as boundary conditions (these were the value 
and slope of the Gibson II potential at r , and the value and slope 
of the truncated Morse potential at r ) - the most general cubic 

(A-3) Ar 3 + Br 2 + Cr + D 

has four unknowns, which just equals the information available. 

The cubic is a standard choice for this "linking potential". 

Solving for the coefficients of the cubic potential 
involved solving four simultaneous linear equations in four un- 
knowns. The technique of solution involved here required putting 
the coefficients of the four linear equations into a matrix which 
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was then transferred to Subroutine CROSYM (developed at Lockheed) 
which performed manipulations on the matrix (using the "Method of 
Crout"). An answer vector was coupled back to the MAIN program. 

The actual potential function used in the program was 
a composite potential, employing: 

the Gibson II potential, for 0 < r < r^, 
the cubic potential, for r^ £ r r 

and the truncated Morse potential, for r < r < r . 

Figure 21 shows the overall composite potential. Note that all 
interatomic separations used in this program were greater than r^, 
that is, the cubic and Gibson II potentials were never used, in 
practice. However, they were included to enchance program flexi- 
bility, since it was anticipated that future improvements of this 
model would allow atomic vibrations which may require the use of 
these other potentials. All these formulas are for interatomic 
spacings in lattice units. A lattice unit is the distance between 
two NNl*s in a single coordinate direction. 

b) . Site Potential Energy. The next step was to assume 
that each lattice site experienced this composite potential as the 
appropriate two-body interatomic interaction potential. Program 
statements 1050 through 1100 set the square of the equilibrium 
interatomic separation with the simple formula 
(A-4) NNDIS2 ( N) = 2*N. 

This formula is exact for the fee, (001) crystal for NNl through 
NN13 in lattice units. Using these values for distance squared, a 
potential was computed for each ith nearest neighbor (i=l,2) based 
on its distance from any particular site. Consider the sketch in 
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Figure 21. The program computed and stored V(r^) and V(r 2 ) . It 
then looked at each site and ( whether it was occupied or not), it 
computed an energy associated with that site based on nearby oc- 
cupied sites out to and including the second nearest neighbor sites. 
This energy was then multiplied by PFIV (a mnemonic for "point 
five" and equal to 0.5 typically, but included as a convenient 
parameter for calculation) since on the average about half the 
potential was associated with each atom of the many two-body pairs 
envisioned. Values of PFIV other than 0.5 gave a different value 
of effective lattice temperature, and in practice lattice tempera- 
ture was varied by changing PFIV. The true equivalent lattice 
temperature based on potential energy only was defined by 

(A-5) EQTEMP = 0 . 5*TEMP/PFIV , 
where EQTEMP is the equivalent lattice temperature 

TEMP is the input temperature used (=1000.0 °K) . 

Figure 22 is a log-log plot of (A-5) . 

In this way each site was given a potential energy 
based on nearby occupied sites. These energies were stored as the 
vector PE(I), where I was the site number. An energy associated 
with both occupied and vacant sites was required in order to deter- 
mine relative transition probabilities from occupied lattice sites 
to nearby vacant sites - and this transition probability was re- 
lated to the energy change between sites by what was essentially a 
Boltzmann factor . 

4 . Transition Probability 

a). Definition. If a given lattice site was occupied and 
had one or more vacant sites as NNl’s, it had (statistically) a 
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finite probability of moving from its present site to the nearby 
vacant site. In practice, when an occupied lattice site was 
reached, an energy difference array was defined between the site 
reached and each of its 12 NNl’s. The array was defined as 
,, DE(J,I) 1 ' - a mnemonic for "difference in energy" - where I was the 
occupied site reached and J (running from 1 through 12) referred to 
the position of one of the NNl’s for site I. The value stored at 
the array position was the energy difference between site I and the 
Jth NN1 position with respect to site I; that is, if K were the 
site number for the Jth NNl position with respect to site I, then 

( A-6) DE(J,I) = PE( K ) - PE(I). 

Since the crystal was always bound, every PE(I) and PE(K) was 
negative, which implied net attraction or crystal cohesiveness. 
However, DE(J,I) may have either positive or negative values, de- 
pending on whether site I or site K was more tightly bound. A 
negative DE(J,I) implied that site K was more tightly bound than 
site I, while a positive DE(J,I) meant that site I was more 
tightly bound than site K. 

Complications arose since it was unphysical to expect 
that any atom (occupied site) had no NNl’s. No atom was permitted 
to jump into a site that was already occupied - this implied destruc- 
tion or loss of matter, which was not allowed. The case of inter- 
changing the atoms between two adjacent occupied sites was 
physically allowable, but uninteresting since it did not change the 
shape of the microstate, and so was not considered. Finally, tran- 
sition from a vacant site was not permitted, since this implied 
creation or addition of matter which was not considered. 
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Sites with no NNl 1 s were not encountered in practice 
(exept for the unphysical "experiment” at 10,000 °K) . The problem 
of jumping from a vacant site and "creating" an atom was easily 
solved by skipping further consideration of a site if it were 
vacant. The cases of jumping from one occupied site to another 
(matter destruction) and of interchanging the atoms in adjacent 
occupied sites (uninteresting interchange) were taken care of by 
defining an arbitrarily large DE(J,I) = 1000.0 in the program 
whenever an occupied site (I) found that the site at its Jth NNl 
position was occupied. Later, when transition probabilities were 
computed, any DE(J,I) ^ 900*0 was ignored, that is, no transitions 
were allowed (or in other words, the transition was given a proba- 
bility of occurrence of 0.0). Sites in the layer z = z 
"looking" out the top of the active lattice volume always "saw" 
site IVAC which was not allowed to be occupied - such sites were 
also given a DE(J,I) = 1000.0 to avoid atom loss. 

Intuitively, physical systems reside at or near the 
bottom of potential wells, so one expects that for 

DE ( J , I ) < 0 

the atom at site I would probably jump to site K, since it would 
then have a lower potential. In fact, by the way transition proba- 
bility was defined, this was exactly what happened. The transition 
probability, p, was defined as [l0]: 

(A-7a) p = 1 DE ( J , I ) ^ 0 

(A -7b) p = exp ( - DE ( J , I ) /kT ) DE(J,I) > 0 

where k = Boltzmann’s constant 

T = absolute temperature. 
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Note that in the program, M kT" was defined at M TAU tT , that is, 

(A-8) TAU = k*T. 

For two or more transitions which satisfy condition 
(A-7a), the computer searched for the deepest well, and carried 
out the transition. Here, that would mean that for DE(J,I) ^ 0 
that the atom at site I jumped to site K. "Carry out the transi- 
tion" meant that the computer zeroed the occupation index of site 
I and placed a "1" in the occupation index of site K. In the 

event that two or more wells had the same (or zero) depth for the 

deepest nearby well, the computer chose one of the sites by a 
random number process. This process involved matrix FTEST which 
was loaded by nesting DO loops at the beginning of the program. 
Matrix FTEST(I,J) was 12 x 12 and had elements 

0.0 below the diagonal, 

1.0 on the diagonal, 

FLOAT ( I ) / FLOAT ( J ) above the diagonal. 

When the energy difference was the same for a jump to two or more 

different sites, the number of sites with the same energy was used 

as J and a DO loop index used to vary I from 1 to J. A random 
number was chosen and the FTEST(I,J) values compared with the 
random number. If the FTEST(I,J) value equalled or exceeded the 
random number, the transition was made to the site which was the 
Ith one which had that particular negative energy difference. If 
the random number exceeded the FTEST ( I, J) value, the loop index 
was incremented and the new FTEST(I,J) value compared with the 
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same random number, and so on, until one FTEST(I,J) value equalled 
or exceeded the random number. In any case, a transition was 
always made for negative (or zero) energy differences. 

The computer decided whether or not two energy dif- 
ferences were "the same" by using a small parameter "EPSLON" - a 
mnemonic for "epsilon". Two energies were considered to be "the 
same" if they were within EPSLON of each other. In practice, 

EPSLON = 1.0 x 10’^. 

For situations in which all DE(J,I) were positive, the 
transition decisions were based on a random number process. First, 
the smallest positive energy difference was found. Then a proba- 
bility for this lowest positive energy difference was computer by 
formula (A-7b) . Next a random number was chosen on the interval 
(0,1). If the transition probability equalled or exceeded the 
random number, the computer carried out the transition, and moved 
on to the next occupied site to repeat the whole process. If the 
probability did not equal or exceed the random number, the com- 
puter found the next lower positive energy difference, calculated 
a probability, and compared it to a new random number. This process 
was continued until either a transition was made to a state of 
higher energy or the site under study ran out of first nearest 
neighbor vacant sites - in this latter case, no transitions were 
made to any other site, and the computer moved on to the next 
occupied site to repeat the procedure for the new site. During 
the calculation, if two or more sites had the same positive energy 
difference, these sites were treated sequentially in a manner simi- 
lar to sites of the same negative energy differences. Positive 
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site energies were considered "the same” if they were within 
EFSLON of each other - this was essentially the same procedure 
used for the negative DE(J,I)'s. 

After any transition, the potentials of the sites, 
affected were "patched". The neighbors of the "old site", that 
is, the previously occupied one, lost an NN1 or an NN2 , therefore, 
PFIV times the proper energy was subtracted from each of these 
NNl f s and NN2 1 s . In addition, the neighbors of the "new site" 
each gained an NNl or an NN2 , so PFIV times the proper energy was 
added to each of these NNl's and NN2 's. 

The question of how appropriate are these probabilities, 
(A-7a) and (A-7b) , may arise. Such probabilities have been used 
in the literature on crystal growth, and have been shown by sta- 
tistical arguments to be correct for an infinitely long chain of 
microstates. According to Leamy and Jackson: 

"averages over a chain of microstates obtained 
Cthrough the use of these transition probabilities] 

... will converge to the equilibrium properties of 
the system in the limit of infinite chain length. "Cio] 

b) . Pseudo-Random Number Generator . The technique used 
to generate random numbers was a standard one called "pseudo- 
random number generation by the multiplicative congruential method". 
The term "psuedo-r andom" implied that the numbers generated were 
not really random, and this was, in fact, true. However, such 
methocfeare generally accepted provided they pass certain statistical 
tests which determine whether the numbers generated are, ibr all 
practical purposes, random. The tests require that the numbers 
generated be uniform on' the interval (0,1), and that the sequence 
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be non-repeating for a sufficiently long period. Other, more 
involved tests also exist. This method has the advantage that it 
is reproducible (so that the same probabilities may be regene- 
rated at will as a check or to show the effect of different initial 
conditions on the final outcome) and fast (this system can work 
"in line") . The term "congruent ial" refers to a congruence re- 
lation in number theory and modulo arithmetic (which deals only 
with remainders) Cl4]. 

Essentially, the method relied on the fact that the 
computer (specifically, the IBM 360/67) had a fixed storage length 
of four bytes (3 2 binary digits or bits) for integer constants. 

This meant that the highest integer the computer could store was 

2 31 - 1 = 2 147 483 647 

a ten digit number (allowing the leading bit for sign plus the 

31 

integer "0", there were 2 -1 locations/combinations for other 

integers [ 3 ]. If one multiplied together two integers, each with 

five or six digits, there was generally an overflow of the machine’s 

storage capacity. FORTRAN did not recognize this as an error, and 

gave only the last eight, nine or ten digits; that is, the leading 

(and generally most significant) digits were truncated, and the 

remainder was stored as an integer. Essentially, it was this 

integer remainder that was the random number. To make it useful, one 

must convert it to a floating point number on the interval (0,1). This 

was done by making the integer a floating point number and dividing 

32 

by 2 , the modulus of the IBM 360 / 67 . Since about half the integer 

remainders (on the average) were negative (due to a negative sign 
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in the sign bit of the integer remainder), one actually had numbers 
on the interval (-0.5, + 0.5), so one must add 0.5 to obtain num- 
bers on the interval (0,1) ['ll]. 

In practice, a total of four FORTRAN statements were 
required to obtain the first random number. They were, for 
example : 

9 KRAN = 16807 

IRAN = 87345 

10 IRAN = IRAN*KRAN 

11 RANDOM = 0.5 + FLOAT ( IRAN ) *2 . 328306 E-10. 

-32 

The number in exponential form in statement 11 is 2 written as 
a decimal number . Whenever one needed a new random number , one 
simply used statements 10 and 11. Note that each random number 
"RANDOM" was generated from the remainder of the previous integer 
remainder ’’IRAN" . This system was fast since it worked ”in line”, 
and computer literature indicates that this generator can produce 
about 7,200 random numbers a second Til]. Since this was a pseudo- 
random generator, its period was finite. However, the theoretical 
maximum period was about 2 or about 537 million numbers Ll4J. 

This type pseudo-random generator is regularly used to generate a 
million or more random numbers. Statistical tests performed on 
the numbers used and observations of the random numbers produced 
indicated that the pseudo-random numbers generated were good enough 
for the practical testing purposes for which they were used here. 

5. Computer Output 

The problem of how to present raicrostate information in 
easily interpreted form was solved by having the computer ’’draw" 
pictures of microcrystal planes, called "microstate pictures”; 
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that is, the computer printed arrays of zeroes and ones to indi- 
cate vacancies or occupied sites, and each digit was in its correct 
x and y coordinate position for the site it represented. See 
Figure 23. The program contained the capability of printing out 
microstate pictures for IX = 10, IX = 20 and IX = 2 6 (with two 
planes side by side), and for IX = any even integer through 60. 

The only restriction on IY and IZ was that they must be even 
integers. The program also contained the capability of choosing 
the appropriate output segment to use. 

Variable numbers of microstates generated were allowed, 
and this was provided via input parameter "NUMRUN" - the number of 
microstates generated. The remainder of the output segment was 
"nice to have" information that was printed out to avoid looking 
through the program. 

For cases in which large numbers of microstates were 
generated, it was desired to have only summary information or 
microstate pictures every mth microstate, and have all microstate 
pictures printed after a certain number of microstates were gene- 
rated. This capability was provided by input parameters Tt MSUM", 
,T MPIX" and ,r MCRIT". MSUM was the index indicating that a short 
summary of the current microstate was printed after each MSUM 
microstates were generated; for example, if MSUM = 5, a summary 
was printed every fifth microstate. The summary was merely the 
number of sites occupied in each z-layer . MPIX was the index 
indicating that a microstate picture was printed after each MPIX 
microstates were generated; for example, if MPIX = 10, a microstate 
picture was printed of every tenth microstate. MCRIT was the index 
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indicating how many microstates were generated prior to printing out 
pictures of every microstate; for example, if MCRIT = 250, every 
microstate from 250 to NUMRUN was printed as a picture. 
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APPENDIX B 



LATTICE LOADING STATEMENTS 

The different initial microstates were generated by sets of 
one or more program statements inserted just prior to statement 
45 in the lattice generator. The statements used for the various 
lattices generated follow. 

1. Perfect Surface 

IF (M-LL/2) 45,45,46 

This statement fills the bottom half of the active lattice 
volume with atoms and leaves the top half empty. 

2. Perfect Surface Plus Randomly Placed Extra Atoms 

IF (M-LL/2) 45,45,49 

49 IF (M-LL/2 - IX*IY/2 ) 48,48,46 

48 IRAN = IRAN * KRAN 

RAND = FLOAT (IRAN) * 2.328306 E-10 
IF (RAND) 45,45,46 

This set of statements fills the bottom half of the active 
lattice with atoms, places atoms randomly on the next plane, and 
leaves all other active lattice planes empty. Sample microstates 
achieved with this set of statements are sketched as Figures 3a 
and 4a . 

3. Perfect Surface Plus Half Plane with Edge Parallel to 

X-Axis 

IF (M - (LL/2 + IX * iy/4) ) 45,45,46 

This statement fills the bottom half of the active lattice 
with atoms , places a half plane of atoms as a monatomic step in 
the next plane, and leaves all other active lattice planes empty. 
The initial microstate achieved is sketched as Figure 5a. 
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4. Pyramids 

IA = 2 
IB = 2 





IF 


(LZ) 20, 


45,20 




20 


IXP 


= LZ* ( IA+1 ) 








IXM 


= IX-1- 


IXP 








IF 


(LX-IXP ) 


46, 


25, 


,25 


25 


IF 


(LX-IXM) 


30, 


30 3 


,46 


30 


IYP 


= LZ* ( IB+1 ) 








IYM 


= I Y-l - 


IYP 








IF 


(LY-IYP ) 


46, 


35; 


,35 


35 


IF 


(LY-IYM) 


45, 


45 j 


,46 



This set of statements loads a pyramid with steps that move 
in three atomic layers at a time - see Figure 6a. To achieve a 
pyramid with steps that move in two (one) atomic layers (layer) at 
a time, change both IA and IB to one (zero) - see Figure 7a. 

5. Pyramid holes 

IA = 2 
IB = 2 

IF (LZ) 20,45,20 
20 IXU = ( IX/2) + LZ* ( IA+1 ) -1 
I XL = (IX/2) - LZ* ( IA+1 ) 

IF (LX-IXL) 45,25,25 
25 IF ( LX-IXU ) 30,30,45 
30 IYU = (IY/2) + LZ*(IB+1)-1 
IYL = (IY/2) - LZ*(IB+1) 

IF (LY-IYL) 45,35,35 
35 IF (LY-IYU) 46,46,45 

This set of statements loads a pyramid hole with steps that 
are three atomic layers wide - see Figure 8a. To achieve a pyramid 
hole with steps two (one) atomic layers (layer) wide, change both 
IA and IB to one (zero) - see Figure 9a. 

6. Ridges 

To achieve a ridge parallel to the y-axis, use the same 
steps as for pyramids, but use the statement 
IB = -1 
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Step width is determined by parameter 1A; for : 



IA = 2 , steps are three atomic layers wide - see Figure 11a; 

IA = 1, steps are two atomic layers wide; 

IA = 0, steps are one atomic layer wide. 

To achieve a ridge parallel to the x-axis , again use the 
same steps as for pyramids, but use the statement 
IA = -1. 

Step width is now determined by parameter IB; for 

IB = 2, steps are three atomic layers wide - see Figure 10a; 

IB = 1, steps are two atomic layers wide; 

IB = 0, steps are one atomic layer wide. 

7. Valleys Parallel to the X-Axis 

IB = 2 

IF (LZ) 30,45,30 
30 I YU = ( I Y/2 ) + LZ* ( IB+1 ) -1 
IYL = (IY/2) - LZ*( IB+1 ) 

IF (LY-IYL) 45 , 35,35 
35 IF (ly-iyu) 46 , 46,45 

This set of statements loads a valley parallel to the x-axis 
with a step width of three atomic layers - see Figure 12 a . To 
achieve such a valley with a step width of two (one) atomic layers 
(layer), change IB to one (zero) - see Figure 13a. 

8. Valleys Parallel to the Y-Axis 

IA = 2 

IF (LZ) 20,45,20 
20 I XU = (IX/2) + LZ* ( IA+1 ) -1 
I XL = (IX/2) - LZ* ( IA+1) 

IF ( LX-IXL) 45,25,25 
25 if (lx-ixu) 46,46,45 

This set of statements loads a valley parallel to the y-axis 
with a step width of three atomic layers '- see Figure l4a . To 
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achieve such a valley with a step width of two (one) atomic layers 
(layer), change IA to one (zero) - see Figure 15a. 

9. (1 -1 1) Plane 

if (lx+lz-ly) 45,45,46 

This statement loads a (1 -1 1) plane - see Figure l6a . 

The plane is truncated if z is smaller than either x or 

^ max max 

y 

max 

10. (-1 1 1) Plane 

IF (LY+LZ-LX) 45,45,46 

This statement loads a (-1 1 1) plane - see Figure 17a. 

This plane is truncated if z is smaller than either x or 
r max max 

y 

max 
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APPENDIX C 



GLOSSARY OF COMPUTER SYMBOLS 

Note: This Appendix also contains the definitions of several terras 



peculiar 


to this simulation. 


A(J,I): 


A 4 x 5 array coupled between the MAIN program and Sub- 
routine CROSYM; it carries information for the computation 
of the coefficients of the cubic potential. 


ALPHA: 


Constant read in as data used in computing parameters for 


BOLTZ : 


the truncated Morse potential. 

The Boltzmann constant, equal to 8.6171 x 10 ^ eV/°K . 



CARRY OUT A TRANSITION: The simulation procedure used to model a 





’’jump”; the occupation index of the original site is set 
equal to zero, and the occupation index of the final site 
(which is an NNl of the original site) is set equal to one. 


CFO: 


A constant computed by the program for use in the force 
equation resulting from the cubic potential. 


CF1 : 


A constant computed by the program for use in the force 
equation resulting from the cubic potential. 


CF2: 


A constant computed by the program for use in the force 
equation resulting from the cubic potential. 


CGBl : 


A constant computed by the program for use in the truncated 
Morse potential. 


CGB2: 


A constant calculated by the program for use in the truncated 
Morse potential. 
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CGDl : 


A constant calculated by the program for use in the trun- 
cated Morse potential. 


CGD2 : 


A constant calculated by the program for use in the trun- 
cated Morse potential. 


CGF1 : 


A constant calculated by the program for use in the force 
equation resulting from the truncated Morse potential. 


CGF2: 


A constant calculated by the program for use in the force 
equation resulting from the truncated Morse potential. 


COMA: 


Labeled COMMON storage. 


CPO: 


A constant computed by Subroutine CROSYM for use in the 
cubic potential. 


CPI : 


A constant computed by Subroutine CROSYM for use in the 
cubic potential. 


CP2 : 


A constant computed by Subroutine CROSYM for use in the 
cubic potential. 


CP3: 


A constant computed by Subroutine CROSYM for use in the 
cubic potential. 


CROSYM : 


A Subroutine developed at Lockheed to solve several 


CVD: 


simultaneous linear equations by the "Method of Crout”. 
"Distance conversion factor", CVR x 10 - ^to convert lattice 
units to meters. 


CVDE : 


CVD/ CVE , a ratio to avoid repeated division; the reciprocal 
of CVED. 


CVE: 


-19 

"Energy conversion factor", 1.6 x 10 to convert elec- 


CVED: 


tronvolts to joules. 

CVE/ CVD, a ratio used to avoid repeated division; the 
reciprocal of CVDE. 
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CVM : 

CVR : 

DCON : 

DE( J,I) : 

DIS : 

DIST( I) : 
DIST2 ( I ) 

EMIN: 



-27 

,T Mass conversion factor", 1.672 x 10 to convert 
atomics mass units to kilograms 

"Radial distance conversion factor", the lattice unit in 
angstroms, used to convert lattice units to angstrom units. 
Constant read/in as data, and used in computing parameters 
for the truncated Morse potential. 

"Difference in energy", an array containing the energy 
difference between an atom at site I and its Jth NN1 
position, whether occupied or vacant (if occupied, the 
particular array value is made arbitrarily large to avoid 
matter destruction); the values on this array are used 
in determining transition probabilities; this array must 
be dimensioned at least 12 x LL. 

SQRT( DIST2 ( I ) ) 

"Distance"; a vector array containing the value of the 
distance (in lattice units) between Ith NN's. 

: "Distance squared"; a vector array containing the value 
of the distance squared (in lattice units squared) between 
Ith NN 1 s . 

"Minimum energy"; the value at the lowest DE(J,I), found 
by searching, when a jump is to be made to a site at lower 
potential; o £ the current value of DE(J,I), found by 
searching and working from lowest to highest energy dif- 
ferences, used when attempting a jump to a site of higher 
potential . 
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EMINXT : "Next minimum energy"; the value of the next lower DE(J,I), 

found by searching, used when attempting to jump to a site 
of higher potential after all attempts at jumping have 
failed for lower positive energy differences. 

EPSLON : "Epsilon"; a constant used to determine if two DE(J,I)’s 

were effectively equal; its value is 1,0 x 10 ^ eV. 
"Equivalent temperature"; the effective lattice temperature 
achieved as a result of varying PFIV; EQTEMP = 0.5 
* TEMP/PFIV , and this relation is plotted in Figure 22. 
"Electron volt". 

A constant read in as data and used in calculation of the 
Gibson II potential. 

A constant read in as data and used in calculation of the 
Gibson II potential. 

Constant read in as data and equal to 2.0*CVR. 

FRCAND(X): "Anderman force”; a function defined in the MAIN program 
for calculation of the force arising due to the truncated 
Morse potential. 

FRCGIB(X): "Gibson force"; a function defined in the MAIN program 
for calculation of the force arising due to the Gibson II 
potential . 

FTEST(I,J): A 12 x 12 array used in choosing the site to which an 
atom jumps if more than one of the NNl sites of an atom 
under consideration have the same energy; for different 
crystal orientations the number of NNl’s may change, and 
therefore the dimensions of array FTEST must also change. 



EQTEMP : 



EV: 

EXA: 



EXB: 



FLAT: 
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FXA: 


A constant computed by the program for use in the force 
equation resulting from the Gibson II potential. 


Is 


A DO loop index. 


IA: 


A constant used by the lattice loader to determine the 
step width of the initial raicrostate for several different 
configurations; see Appendix B. 


IB: 


A constant used by the lattice loader to determine the 
step width of the initial microstate for several different 
configurations; see Appendix B. 


I COUNT : 


A counting index which is the number of the raicrostate 
being generated by the program. » 


IFUL : 


"Full"; the number of a lattice site in the perfect cry- 
stalline substrate on which the active lattice volume is 
placed; all lattice sites in this substrate are always 
full and each is represented by site number IFUL* 


IH2 : 


Alphanumeric array read in as data for truncated Morse 
potential parameters. 


IHEAD: 


Alphanumeric array read in as data to provide an output 
heading for the particular simulation attempted. 


IHNOl : 


Alphanumeric array read in as data to provide a program 
flag warning that during a calculation an atom was found 
that had no NN1 f s (this is unphysical in this simulation); 
the message merely states "NO NNl ? s". 


IHTARG : 


Alphanumeric array read in as data and used to define the 
particular metal simulated. 


IHTPOT : 


Alphanumeric array read in as data to provide a title for 
Gibson II potential parameters. 
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II: 

IJ: 

IK: 

IRAN: 

IRANI : 

ISX: 
ISXM : 

ISXP : 

ISZ: 

ISZM : 

ISZP: 

IT: 

ITEST: 



ITT: 



A DO loop index. 

A DO loop index. 

A DO loop index. 

The present (or initial) value of the constant used in 
computing pseudo-random numbers. 

"Initial IRAN"; defined to retain a knowledge of the 
original random number seed used in a simulation. 

IX/2, a shift index indicating a shift of one row of atoms. 
ISX-1, a shift index indicating a shift of one atom less 
than a full row. 

ISX+1, a shift index indicating a shift of one atom more 
than a full row. 

ISX*IY, a shift index indicating a shift of one plane of 
atoms . 

ISZ-1 , a shift index indicating a shift of one atom less 
than a full plane. 

ISZ+1 , a shift index indicating a shift of one atom more 
than a full plane. 

An index used in computing ITT. 

A IX) loop index*, also, a constant used to determine which 
site may receive an atom when carrying out a transition 
in the case of several available sites with the same 
energy . 

IT+JT+KT, the index which determines when or not a particular 
coordinate position should be labeled as a valid lattice 
site by the lattice generator. 
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I VAC: 



IX: 



IXL : 



IXM : 



IXP: 



IXU : 



IY: 



IYL: 



IYM : 



IYP : 



I YU: 



IZ: 



J: 



’’Vacant"; the number of a lattice site in the void above 
the active lattice volume; all sites in this void are 
always vacant and each is represented by site number IVAC. 
An even integer read in as data indicating the number of 
active lattice volume planes generated in the x-direction. 
A variable used as a shifter when generating pyramid holes 
and valleys; see Appendix B. 

A variable used as a shifter when generating pyramids and 
ridges; see Appendix B. 

A variable used as a shifter when generating pyramids and 
ridges; see Appendix B. 

A variable used as a shifter when generating pyramid holes 
and valleys; see Appendix B. 

An even integer read in as data indicating the number of 
active volume planes generated in the y-direction. 

A variable used as a shifter when gener at ing . pyramid holes 
and valleys; see Appendix B. 

A variable used as a shifter when generating pyramids and 
ridges; see Appendix B. 

A variable used as a shifter when generating pyramids and 
ridges; see Appendix B. 

A variable used as a shifter when generating pyramid holes 
and valleys; see Appendix B. 

An even integer read in as data indicating the number of 
active lattice volume planes generated in the z-direction. 
A DO loop index. 
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JJ: 


IY-J, a constant used in the output segment to draw and 
label microstate pictures. 


JJM : 


JJ-1, a constant used in the output segment to draw and 
label microstate pictures. 


JLOE(I) : 


"Low energy index"; a vector array used to record the site 
number of the lowest energy NN1 site of an atom attempting 
to make a transition; this vector must have dimension of 
at least 12. 


JMAX: 


J-l, an index indicating the number of NN1 sites of lower 
energy for an atom attempting to make a transition. 


JP1 : 


J+l, a constant used in the output segment to draw and 
label raicrostate pictures. 


JT: 


An index used in computing ITT. 


JUMP: 


The physical process of an atom leaving its present 
location and moving into a vacant NNl site. 


KA: 


K*ISZ, a constant used in the output segment to draw and 
label microstate pictures. 


KB: 


KA-J*ISZ+1, a constant used in the output segment to draw 
and label microstate pictures. 


KC: 


KA-JP1*ISX+1 , a constant used in the output segment to draw 
and label microstate pictures. 


KD: 


KB+ISXM, a constant used in the output segment to draw and 
label raicrostate pictures. 


KE: 


KC+ISXM, a constant used in the output segment to draw and 
label raicrostate pictures. 


KF: 


KP1*ISZ, a constant used in the output segment to draw and 
label microstate pictures. 
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KG: 


KF-J*ISX+1, a constant used in the output segment to draw 
and label microstate pictures. 


KH: 


KF-JP1*ISX+1 , a constant used in the output segment to 
draw and label microstate pictures. 


KI: 


KG+ISXM, a constant used in the output segment to draw 
and label microstate pictures. 


KJ: 


KH+ISXM, a constant used in the output segment to draw 
and label microstate pictures. 


KKOUNT : 


An index used to avoid recalculating potential coefficients 
when more than one crystal simulation is attempted. 


KM1 : 


K-l,an index used in printing out raicrostate summary infor- 
mation; also, a constant used to draw and label raicrostate 
pictures. 


KNEW: 


"New K"; site to which a transition has just been carried 
out . 


KOUNTM : 


K0UNT2-1, an index used in deciding which site to jump 
to when there is more than one NN1 site of lowest energy 
to choose from. 


KOUNT1 : 


A counting index used to determine the number of NN1 sites 
that have lower or higher energy with respect to the site 
under consideration. 


KOUNT2 : 


A counting index used to determine the multiplicity of 
NN1 sites with lowest energy. 


KOUNT 3 : 


A counting index used to determine the multiplicity of 
NN1 sites with higher energy. 


KOUNT4: 


A counting index used to determine whether or not all higher 
energy sites have been tested for the possibility of making 
a transition. 
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KOUNT 5: 


An index controlling program flow when attempting to 
determine if all sites of higher energy have been tested. 


K0UNT6 : 


A counting index indicating the multiplicity of NN1 sites 
with next higher energy. s 


KP1 : 


K+l , an index used in printing out raicrostate summary 
information; also, a constant used to draw and label 
microstate pictures. 


KP2: 


K+2, an index used in printing out microstate summary 
inf or mation . 


KP3: 


K+3, an index used in printing out microstate summary 
inf or mat ion . 


KRAN: 


A constant used in computing pseudo-random numbers. 


KT: 


An index used in computing ITT. 



LATTICE UNIT(L.U.): The distance between two NN1 1 s along a single 





coordinate direction. 


LDX: 


A shift index used to assign lattice site coordinate 
positions in the x-direction. 


LDY : 


A shift index used to assign lattice site coordinate 
positions in the y-direction. 


LDZ: 


A shift index used to assign lattice site coordinate 
positions in the z-direction. 


LL: 


"Lattice length"; the number of available sites in the 
active lattice volume, equal to IX*IY*IZ/2. 


L.U : 


See LATTICE UNIT. 


LX: 


The x-coordinate position assigned to a particular site 
during lattice generation. 
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LY: The y-coordinate position assigned to a particular site 

during lattice generation. 

LZ: The z-coordinate position assigned to a particular site 

during lattice generation. 

M: A counting index used to assign a label to lattice sites; 

also, an index (= MCRO) used in Subroutine CROSYM. 

MCRIT: A constant read in as data; all microstates generated from 

microstate number MCRIT to NUMRUN are printed out as 
microstate pictures. 

MCRO: A constant read in as data, indicating the number of 

unknowns to be solved for by Subroutine CROSYM. 

MICROSTATE: The surface configuration either initialized in the 
program or generated by it. 

MICROSTATE PICTURE: An array of occupation indices printed out in 
their correct coordinate positions to give a physical 
picture of the surface configuration (see Figure 23) . 

MICROSTATE SUMMARY: A listing of the microstate number and the 
number of atoms in each z -plane. 

MPIX: A constant read in as data indicating that between micro- 

state zero and microstate MCRIT, a raicrostate picture is 
printed out every MPIXth microstate. 

MSUM: A constant read in as data, indicating that between micro- 

state zero and microstate MCRIT, a microstate summary is 
printed every MSUMth microstate (except when a microstate 
picture is printed). 
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MULT(I): Multiplicity"; a vector array indicating the numbers of 



NNl sites which have the same energy difference with 
respect to the atom of interest; this vector must have 
1 dimension at least as great as the number of NNl's for a 
lattice location (in this simulation, the proper number 
was 12) . 

NA: An index indicating which site anatom jumped to when 

several NNl sites of the same energy were located; the 
value of NA was chosen stochastically. 

NBRFOR ( J ,K ) : "Neighbor four"; an array containing the site numbers 
of the 12 NNVs of the Ith lattice site; since this simu- 
lation did not use NN4 f s the array was dimensioned lxl, 
but if NN4 's are included, it must have dimensions of at 
least 12 x LL. 

NBRONE(J,I): "Neighbor one"; an array containing the site numbers 
of the 12 NNl f s of the Ith lattice site; this array must 
have dimensions of at least 12 x LL. 

NBRTHR(J,I): "Neighbor three"; an array containing the site numbers 
of the 24 NN3 1 s of the Ith lattice site; since this simu- 
lation did not use NN3's, the array was dimensioned lxl, 
but if NN3*s are included, it must have dimensions of at 
least 24 x LL. 

NBRTWO(J,I): "Neighbor two"; an array containing the site numbers 
of the six NN2 f s of the Ith lattice site; this array must 
have dimensions of at least 6 x LL . 
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NCRYST: An index read in as data to indicate the number of crystal 

simulations to be run using a single input deck; if 
NCRYST = 0, the simulation with which this value is read 
in is the last (or the only) simulation run; if NCRYST j* 0, 
there is at least one simulation following the one con- 
taining a non-zero value of NCRYST. 



NN: 


"Nearest 


neighbor " . 


NNl: 


"First nearest neighbor". 


NN2: 


"Second 


nearest neighbor". 



NNDIS2(I): "Nearest neighbor distance squared"; the square of the 

distance to the Ith NN; this vector must have dimension 
at least as great at NUMNN. 

NNFOOC(I): "Number of NN4 sites occupied", defined with respect to 
site I; this vector must have dimension at least as great 
as LL. 

NNONOC(I): "Number of NNl sites occupied", defined with respect 
to site I; this vector must have dimension at least as 
great as LL. 

NNTHOC(I): "Number of NN3 sites occupied", defined with respect to 

site I; this vector must have dimension at least as great 
as LL. 

NNTWOC(I): "Number of NN2 sites occupied", defined with respect to 
site I; this vector must have dimension at least as great 
as LL. 

NOCC(I): Vector array containing the occupation index of each 

lattice site; this vector must have dimension at least 
as great as LL+3 . 
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NPLANE: "Plane number”; NZ(I)+1, a constant used to define plane 
numbers greater than zero for use as subscripts. 

NUMATM : “Number of atoms”; the number of atoms actually placed in 

the active lattice volume by the lattice generator. 

NUMNN : "Number of nearest neighbors”, considered in a particular 

simulation . 

NUMRUN: "Number of runs”; a constant read in as data indicating 

the number of microstates to be generated by the simu- 
lation . 

NUMSIT(I); "Number of sites" occupied in the Ith plane; here, I is 
defined by NPLANE in order to have non-zero subscripts. 

NX(I): Vector array containing the values of the x-coordinate 

for each lattice site; this vector must have dimension 
at least as great as LL. 

NY(I): Vector array containing the values of the y-coordinate 

for each lattice site; this vector must have dimension 
at least as great as LL. 

NZ(I): Vector array containing the values of the z-coordinate 

for each lattice site; this vector must have dimension 
at least as great as LL. 

OCCUPATION INDEX: A number indicating whether or not a site is 
occupied; an occupation index of zero implies that the 
site is vacant, and an occupation index of one implies 
that the site is occupied; no other numbers are allowed 
as an occupation index. 



73 



PE(I): "Potential energy”; a vector array containing the value of 



the potential energy at each lattice site due to nearby 
occupied sites; this vector must have dimension at least 
as great as LL . 



PENN(I) : 


"Potential energy due to nearest neighbor"; the potential 
energy contribution due to the Ith NN. 


PFIV : 


"Point five"; an energy distribution factor employed in 
calculating site potential energy and equivalent lattice 
temperature; the correct average value of PFIV is 0 . 5 . 


POTF(X) : 


A function defined in the MAIN program for calculation 
of the Gibson II potential. 


P 0 T 2 F( X) 


: A function defined in the MAIN program for calculation 
of the cubic potential. 


P 0 T 3 F(X) 


: A function defined in the MAIN program for calculation of 
the truncated Morse potential. 


RAND: 


FLOAT ( IRAN) * 2 . 328306E-IO , a pseudo-random number on the 
interval (-0.5, 0.5) used in placing atoms randomly on a 
perfect surface; see Appendix B. 


RANDOM: 


0 . 5 +FL 0 AT( IRAN) * 2 . 328306E-IO , a pseudo-random number on 
the interval (0,1) used for comparison with probabilities 
when a decision based on a probability was made. 


RE: 


A constant read in as data used in computing truncated 
Morse potential parameters. 


ROEA: 


A constant read in as data equal to the distance in 
lattice units at which the Gibson II and cubic potentials 
match . 
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ROEB: 

ROEC: 

ROEC2 : 

SUM: 

TAR: 

TAU: 

TEMP: 

TMAS : 

TPOT : 

TPROB : 



A constant read in as data equal to the distance in 
lattice units at which the cubic and truncated Morse 
potentials match. 

A constant read in as data equal to the distance in 
lattice units at which the truncated Morse potential is 
truncated. 

"ROEC squared”; ROEC*ROEC. 

A variable used by Subroutine CROSYM . 

Alphanumeric array read in as data and not used in this 
s imulation . 

BOLTZ*TEMP, i.e, the ”kT" factor used in computing a 
Boltzmann factor. 

"Temperature"; the nominal lattice temperature; see 
EQTEMP . 

"Mass"; a constant read in as data equal to the mass of 
a lattice atom in atomic mass units. 

"Total potential"; the sum of all occupied site potentials 
due to interactions with nearest neighbors through NUMNN. 
"Transition probability", defined by equations A-7. 
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COPPER CRYSTAL SURFACE DYNAMICS SIMULATION 



COMPUTER PROGRAM 
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GO TO 195 

SPECIAL CASE: X=IX-1, Y=0, Z=IZ-1, TOP RIGHT FRONT CORNER 
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TABLE I 



POTENTIAL FUNCTIONS AND COEFFICIENTS 



V G ibs on ll R > = exp(EXA + EXB*R) 
v cubic ( R ) = R*( R *( R * cp 3 + CP2)+CPl)+CPO 

V trun Morse (R) = exp(CGDl+CGBl*R) - exp(CGD2 +C GB2*R) 



Type Potential 


Coefficients 


Cutoff Distance 
( latt ice units ) 


Gibson II 


EXA = 10.02407 
EXB = -9.19667 


r A = °' 83 


cubic 


cpo = 587.6182 
CPI = - 1593.863 

CP 2 = 1450.286 

cp 3 = - 442.2820 


1 B = la ° 


truncated Morse 


CGD1 = 6.649136 
CGD2 = 3.651313 
CGB1 = -5.077808 
CGB2 = -2.538904 


r c = 2.40 
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RESULTS OF A SINGLE ATOM JUMPING ON A PERFECT SURFACE 
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Vertical Motion Only Lateral Motion Only No Motion 




Figure 1. The Active Lattice Volume. 




Figure 2. Overall Program Flow. 
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Figure 3 . Perfect Surface Plus 99 Randomly Placed Atoms. 

(a) Initial raicrostate (lines bordering atoms 
and atom clusters are parallel and perpendicular 
to the { lOO) direction). 
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Figure 4. Perfect Surface Plus 104 Randomly Placed Atoms. 

(a) Initial microstate (lines bordering atoms 
and atom clusters are parallel and perpendicular 
to the { lOO) direction). 
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Figure 5. Perfect Surface Plus Half Plane Monatomic Step; 
(a) initial microstate 




Figure 5b. Microstate 500 (lines bordering the atom 
cluster are parallel and perpendicular 
to the (lio) direction) 
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Figure 6. 



Pyramid With Steps Three Atomic Layers Wide; 
(a) initial raicrostate. 




Figure 6b. Microstate 500 (lines bordering the atom 
cluster are parallel and perpendicular to 
the (lio) direction). 
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Figure 7 . Pyramid With Steps Two Atomic Layers Wide, 
(a) Initial raicrostate. 




(b) Microstate 500 (lines bordering atom clusters are 
parallel and perpendicular to the (110) direction); the 
hole in the z - 1 plane is a result of the partial fill-in 
and reorientation of the initial grooves in the z = 1 plane 
however, the vacancy in the z = 0 plane (lower left corner) 
is due to an atom jumping out of a perfect surface (the 
z = 0 plane) - this phenomenon was observed infrequently. 
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Figure 8. Pyramid Hole With Steps Three Atomic Layers 
Wide; (a) initial microstate. 




Figure 8b. Microstate 500 (lines bordering atoms and 

atom clusters are parallel and perpendicular 
to the (110) direction). 
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Figure 9. Pyramid Hole With Steps Two Atomic Layers Wide, 
(a) Initial microstate. 




(b) Microstate 500 (lines bordering clusters of atoms are 
parallel and perpendicular to the (lio) direction); the 
holes appearing here are due to the partial fill-in and 
reorientation of the initial holes. 
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Figure 10. Ridge Parallel to X-Axis With Steps Three 

Atomic Layers Wide, (a) Initial microstate. 




(b) Microstate 500 (lines bordering atom clusters are 
parallel and perpendicular to the <110> direction); the 
hole appearing in the z = 1 plane is due to the partial 
fill-in and reorientation of the initial groove in the 
z = 1 plane. 
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Figure 11. Ridge Parallel to Y-Axis With Steps Three 

Atomic Layers Wide, (a) Initial microstate. 




(b) Microstate 500 (lines bordering atom clusters are 
parallel and perpendicular to the (lio) direction); the 
hole appearing in the z = 1 plane is a result of the 
partial fill-in and reorientation of the initial groove 
in the z = 1 plane. 





(c) Microstate 451 , using a different random number seed 
then used to achieve the microstate depicted in Figure lib 
(the lines bordering atoms and atom clusters are parallel 
and perpendicular to the ( lio) direction); the hole 
appearing in the z = 1 plane is a result of the partial 
fill-in and reorientation of the initial groove in the 
z = 1 plane. 
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Figure 12. Valley Parallel to X-Axis With Steps Three 
Atomic Layers Wide; (a) initial microstate. 




Figure 12b. Microstate 500 (lines bordering the atom 

clusters are parallel and perpendicular to 
the <110> direction). The hole appearing 
here is a result of the partial fill-in and 
reorientation of the initial groove in the 
2=1 plane. 
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Figure 13 . Valley Parallel to X-Axis With Steps Two 

Atomic Layers Wide; (a) initial microstate. 




Figure 13b. Microstate 500 (lines bordering the atom 

clusters are parallel and perpendicular to 
the (110> direction). The hole appearing 
here is a result of the partial fill-in 
and reorientation of the initial groove in 
the 2,-2 plane. 
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Figure 14. Valley Parallel to Y-Axis With Steps Three 
Atomic Layers Wide, (a) Initial microstate. 




(b) Microstate 500 (lines bordering atom clusters are 
parallel and perpendicular to the ( 110) direction); the 
hole appearing in the z = 1 plane is a result of the 
partial fill-in and reorientation of the initial groove 
in the z = 1 plane. 




Figure 15. 



Valley Parallel to Y-Axis With Steps Two Atomic 
Layers Wide, (a) Initial raicrostate. 




(b) Microstate 500 (lines bordering atom clusters are 
parallel and perpendicular to the <110> direction). 
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Figure 17. (-1 1 1) Plane, Truncated, (a) Initial 
microstate. 




(b) Microstate 500. 
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(C) 
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(b) 



Figure 18. The 12 NN1 Positions. 

(a) In the z-plane above 
the atom; 

(b) In the same z-plane as 
the atom; 

(c) In the z-plane below 
the atom. 
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Figure 19. The Six NN2 
Pos itions . 

Sites at positions 1-4 are 
in the same z-plane as the 
atom; position 5 is two 
z-planes below the atom; 
position 6 is two z-planes 
above the atom. 
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EQTEMP 

(°K) 




Figure 22. Equivalent Lattice Temperature as a 
Function of PFIV. 
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Z = 0 PLANE 

NUMBER OF SITES OCCUPIED = 10 



Figure 23. Sample Microstate Picture. Here, IX = 10 
and IY = 4. 
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Figure 24. Microstate 310 of an Initially Perfect Surface 
at 10,000 °K. (a) Lines bordering atoms and 

clusters of atoms are parallel and perpendicular 
to the (lOO) direction. 
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Figure 25. Microstate 180 of an Initially Perfect Surface 
at 2,000 °K (lines bordering atom clusters are 
parallel and perpendicular to the (lio) 
direction). Note that atoms have jumped out of 
the perfect surface and coalesced on the surface. 




Figure 26. Microstate 190 of an Initially Perfect Surface 
at 2,584 °K (lines bordering atom clusters are 
parallel and perpendicular to the (lio) 
direction). Note that a hole two atomic layers 
deep has been burroughed as atoms left the 
perfect surface and piled up near it. 
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Figure 28. Break Repair Sequence For Atom Aggregates in a 
Single Plane. The initial microstate was Figure 
10a. 
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Figure 28f. Microstate 150, z = 2 plane (63 atoms). The 
two groups of atoms are held together only 
by one NN2 interaction. 
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Figure 28g. Microstate 170, z = 2 plane (6l atoms). 





Figure 28i. Microstate 190, z = 2 plane (60 atoms). The 
two groups of atoms are joined only by one 
NN2 interaction. 
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Figure 2 8 j. Microstate 200, z = 2 plane (59 atoms). 





x 



Figure 28k. Microstate 2ZfO, z = 2 plane (57 atoms). 




Figure 28l. Microstate 290, z = 2 plane (63 atoms - i.e., 
some atoms have fallen down from the z = 3 
plane) . 
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Figure 28m. Microstate 330, 2 = 2 plane (6l Atoms). 





Figure 28n. Microstate 350, z = 2 plane (57 atoms). 





Figure 280. Microstate 370, z = 2 plane (55 atoms). 



Y Y 




(a) z even (b) z odd 

Figure 29. Relative Atomic Positions in a Lattice Plane. 
Here, IX = 4 and IY = 6. 
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(a) z = 0 



(b) z = 1 



Figure 30. 



The Site Numbering Scheme Imposed by the 
Lattice Generator. Here, IX = 4 and IY = 6 . 
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(a) Overview. 
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(b) Specific site, images for IX = IY = 6 , 
z even. 

Figure 31. Visualization of Periodic Boundary Conditions. 
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